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^ ; Abstract 

. ■ In this paper, the imaginary-time path integral representation of the canonical partition function 

Cd ' of a quantum system and non-equilibrium work fluctuation relations are combined to yield methods 

^ . for computing free energy differences in quantum systems using non-equilibrium processes. The 

i-rt ! path integral representation is isomorphic to the configurational partition function of a classical 

^ \ field theory, to which a natural but fictitious Hamiltonian dynamics is associated. It is shown that 

if this system is prepared in an equilibrium state, after which a control parameter in the fictitious 

Hamiltonian is changed in a finite time, then formally the Jarzynski non-equilibrium work relation 

^S| ■ and the Crooks fluctuation relation are shown to hold, where work is defined as the change in 

'^ \ the energy as given by the fictitious Hamiltonian. Since the energy diverges for the classical field 

f~^ \ theory in canonical equilibrium, two regularization methods are introduced which limit the number 

'nI" ■ of degrees of freedom to be finite. The numerical applicability of the methods is demonstrated for 

• ' a quartic double-well potential with varying asymmetry. A general parameter-free smoothing 

^—v . procedure for the work distribution functions is useful in this context. 

OO ■ 
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I. INTRODUCTION 

A number of relations valid in the far from equilibrium regime have appeared in the last 
fifteen years [l|, y, yl, U, |5|, la, [Zl, |8|] that show intriguing relationships between fluctuations in 
non-equilibrium systems governed either by deterministic or stochastic dynamics. Among 
these relations, the JarzynskifsJ, |6|] and Crooks relations [7|, |8| provide a means to compute 
the free energy difference between two classical systems by use of a control parameter that 
switches the system from one ensemble to another in a well-defined manner. The extension of 
these relations to quantum systems has been analyzed recently by several authorsJ9l.llCll.llll. 
12|, [l3|, [ij] . Current methods of calculating free energy differences in quantum systems using 
non-equilibrium processes rely on the knowledge of the quantum history of the system which, 
while yielding a conceptually appealing picture, cannot provide a reasonable scheme for the 
computation of free energy differences in practical applications. The challenges associated 
with constructing the correct coherent quantum dynamics make the approach difficult to 
implement. 

The path integral representation of the canonical partition function is based on mapping 
a quantum system at finite temperature onto a classical system with additional degrees of 
freedom [l5|,[l6|,[l7|. A non-equilibrium process can be carried out on this isomorphic classical 



system along a well defined trajectory in fictitious time. As will be demonstrated, the 
Jarzynski and Crooks relations are valid for such a process provided conditions analogous to 
the ones required for the relations in a classical system are met. As a consequence, quantum 
free energies using fictitious non-equilibrium classical processes can be obtained using the 
path integral representation. The path integral formulation, however, involves an infinite 
number of degrees of freedom, and, as a result, non-equilibrium dynamical processes in this 
representation could lead to divergences which need to be regularized. The practicality of 
the method depends sensitively on the rate of convergence of properties determined from the 
regularized path integral to their true quantum values, an issue that is given considerable 
attention in this paper and the following paper in the series, which treats the case of a 



harmonic oscillator in detail 18 



The paper is organized as follows. In Sec. [ITl the path integral formalism and two 
equilibrium methods to compute free energy differences are reviewed. The Jarzynski non- 
equilibrium work relation and the Crooks fluctuation relation using fictitious non-equilibrium 
processes are derived in Sec. IIIII In Sec. HVl the subtleties associated with the divergence 
of the total energy and work of a system with an infinite number of degrees of freedom 
are discussed, and two regularization methods are introduced. The first one is based on 
the Fourier representation of closed paths representing quantum particles and the second is 
based on a spatial discretization of the closed paths. In Sec. |Vl the non-equilibrium free 
energy method is applied numerically to the free energy of a double well potential of quartic 
form as the asymmetry between the wells is varied. The conclusions are given in Sec. I VII 

II. SYSTEM AND DEFINITIONS 

Consider a quantum system with a Hamiltonian of the form H{X) = T + V with T = 
fp'/{2m) and V = V{x,X). For simphcity of presentation, the position operator x and the 
associated momentum operator p here are taken to be one-dimensional though the extension 
of the analysis to higher dimensional systems is straightforward. Note that the potential 
energy V depends on a control parameter A which is independent of the configuration of 



the system. The canonical partition function of this system at an inverse temperature (3 is 
defined by 

Z(A) = Tre-^^(^), (1) 

which is related to the free energy by 

Z(A)=e-^^W. (2) 

The partition function can be written in path integral form as 15|, [la 



Z{X)= fvxe-^^^''^\ (3) 

where the integral is over closed paths x{s) [i.e., x{[3h) = x(0)] and the Euclidean action S 
is a functional of x given by 



S[x, A] = / ds 
Jo 



(4) 



1 / dx \ *^ 



Here and below the s dependence of x in integrals over s will always be implied. Due 
to the cyclic property of the trace, the quantum mechanical equilibrium ensemble average 
(A)^ = TT{Aexp[—(3H{X)]}/Z{X) of an operator A{x) can be written in the path integral 
formulation as 

{A)T = ^ jvxAlx, A] e-^^[^'^] = (A[x, A]),, (5) 

where (. . .)x denotes a path integral average -^kr jVx . . . e~E'^[^'''^l and A denotes the imag- 
inary time average 

A[x,X] = ^J dsA{x{s),X). (6) 

Note that the path integral average of a function of a single imaginary time point, 
{A{x{s*))) X, does not depend on the choice of the time point s*, due to the imaginary time 
translational invariance of the Euclidean action. One may therefore replace {A{x{s*))) x by 
the imaginary time average (^[x]) a, which is often advantageous for reasons of computational 
efficiency. 

Eq. ([5]) defines the equilibrium ensemble average of a general functional A[x, X] of the path 
x{s), which could also be in general a function of A. Note that for multiple-particle systems, 
one needs to incorporate exchange effects associated with quantum statistics |l5|, [ly . 



In the context of equilibrium statistical mechanics, the free energy difference between 
two systems characterized by different values of the control parameter A can be computed 
in this picture through the so-called thermodynamic integration method, and is given by 

7a, \pnj, dx 

^ .a(^|..A|)^. (7) 

where the free energy change refers to the difference between the initial and final states 
through A (i.e. AF = F{Xb) — F{Xa))- Alternatively, one can compute the free energy 



change directly using the following identities [15 



AF = — — log(exp 

P 

1 



--{S[x,Xb] - S[x,Xa]) 
n 



j log (exp [-/5 ( V[x, Xb] - V[x, Xa])] >,^ . (8) 

In the classical limit, the closed paths x{s) transform into a point x and Eqs.Q and (IE]) 
transform into the well-known classical expressions [l9l. |20 . 



III. NON-EQUILIBRIUM RELATIONS 

The non-equilibrium relations derived in this section pertain to a non-equilibrium process 
in fictitious time. To construct the fictitious dynamics, a new field p{s) is introduced which 
is also periodic in imaginary time, satisfying p{s) = p{s + j3h). By multiplying the path 
integral representation of the partition function in Eq. ([3]) by 



Vpe s/o^^d^li 



C" 



(9) 



where C is a normalization constant and /i is an arbitrary fictitious mass, one obtains 



Z{X) = C VxVpe 



-i5M]-i/d. 



2/j 



(10) 



Using Eq. (j4]), this equation can be cast in the form of a classical partition function 

Z(A) = C [vxVpe-^"^''^P^^\ (11) 



where the fictitious Hamiltonian is given by 



H[x,p,X] 



du 



p^ 
2/1 



^<l)' + "-(^-^' 



(12) 



Here, a scaled imaginary time variable u = s/{Ph), has been introduced, while 

m 

Equations (ITTl) and flT^ correspond to the thermal field theory of a classical closed elastic 
string of unit length in one dimension with mass /i and string tension k. 

Recently, Scholl-Paschinger and Dellago showed that Jarzynski's non-equilibrium work 
relation [sl, |6| holds for a wide class of classical deterministic systems with a finite number 
of degrees of freedom 21j . It was demonstrated that if the dynamics of a system for fixed 
values of A admits an invariant distribution of the system plus bath equal to the canonical 
distribution at an inverse temperature /? multiplied by a function dependent on the bath 
variables only, then the free energy difference can be found from 



-/3AF 



(exp(-/3iy))A„ 



(13) 



where the average is over repetitions of a non-equihbrium process in which the system starts 
from a configuration drawn from the above mentioned invariant distribution at A = A^ and 
is driven out of equihbrium by varying the control parameter A in a finite amount of time r 
from A^ to A^ via an arbitrary protocol A(t). Furthermore, the work W done in the process 
in Eq. fll3p is given by 

The result of Scholl-Passinger and Dellago holds for a variety of different deterministic 
dynamics. Thus, for the purpose of computing AF from Eq. (fT3|l . any dynamical evolution 
scheme can be used. 

Perhaps the simplest dynamical evolution is generated by Hamiltonian dynamics. In the 
current context, this dynamics is governed by the fictitious Hamiltonian in Eq. (fT2|) . so the 
work W is also fictitious. The equations of motion for the fields x{u, t) and p{u^ t) resulting 
from the fictitious Hamiltonian are 

dx ^ 5H[xpX] ^ v_ (,4^) 

dp 5H[x,p,\] d'^x d 

m = -^^(^ = '^a^ - a^^^"' ^)- ^^'^^ 

Equations (I14al) and fll4bl) are the usual equations of motion of a single elastic string in an 
external potential V. It should be stressed that these equations have no relation to the real 
time evolution of the original quantum particle. 

Because the system is isolated, the fictitious work done on the system in changing the 
control parameter A is precisely the difference between values of the fictitious Hamiltonian 
at time r and at time 0. Introducing the convention that quantities without explicit time 
arguments are taken at time zero, one can thus write 

W = H[x{t),p{t),Xb] - H[x,p,Xa]. (15) 

Now consider the exponential average of W over initial conditions drawn from the canonical 
equilibrium of the string at A = Ayi: 

,pw\ _ JVxVpe-f^^e-l^^^^'P'^^^ 
fVxVp e~^^[^W'PW'^s] 



JVxVp e-^^i^^'P'^A] 



In the numerator, one can change path integration variables from the initial field x and p 
to x' = x{t) and p' = p{t). The Jacobian of this transformation is equal to unity due to 



Liouville's theorem |22[|. so that 

JVx'Vp' e-/^^[^''P''^s] 



(e-^^> 



e-^^^, (16) 



Z{Xb) _ -f3AF 



Z{K 



which is Jarzynski's non-equilibrium work relation, i.e. Eq. flT3|) . It should be stressed that 
while W is fictitious work, the resulting AF is the real quantum free energy difference. 



Note that for a process that occurs infinitely fast, i.e. the switching time r = 0, one 
recovers Eq. ([8]). For an infinitely slow process, however, one doesn't recover the canonical 
equilibrium expression Eq. (jTj) due to the fact that the system evolves in isolation. 

Another non-equilibrium relation, the Crooks fiuctuation relation[7|, l8|, can also be shown 
to hold in this context. Consider a non-equilibrium process performed as stated above, as well 
as in the reversed sense, i.e. starting from configurations drawn from a canonical distribution 
at A = As and driven out of equilibrium by varying A in time from A^ to A^ in the reversed 
direction of time, that is, with A(t) —^ A(r — t). Then the probability that an amount W of 
work is done during the forward process can be written as 

PfiW) 

-I3H[x,p,Xa] / r-r _ f)ZJ 

VxVp 5 iW- dt X—- 

Za V Jo 9X 

Vx'Vp' 6 iw+ dt A^ 

Za V Jo 9X 

= ef'^e-^^^Pri-W), (17) 

where Pr{—W) is the probability that an amount of work —W is done during the reverse 
process. Eq. ( IT71) is known as the Crooks fiuctuation relation|7|, \^. Note that the value of 
W at which the two distributions Pf{W) and Pr{—W) become equal, which will be denoted 
by Wc, is precisely when Wc = AF. Thus this relation allows AF to be computed by 
determining where the plots of Pf(W) and Pr{—W) versus W intersect. This approach is 
known as the crossing method[23|, |2J]. 

The Crooks fiuctuation relation, Eq. fITTj) . can be extended to a conditional ensemble in 
which a particular value of a variable, called the "reaction coordinate" , is held fixed to a 



value x[25[. The free energy f{x) at this constrained value of the reaction coordinate is 
known as the potential of mean force. Here, we consider the reaction coordinate x to be 
given by a functional x of the position- field x, and the Hamiltonian to take the form 

H[x,p,X] = Ho[x,p] + ^xHX). (18) 

The potential of mean force f{x) is then defined as 

^-mx) = (^s{x-xM))o (19) 

plus an arbitrary constant, where the average is over the path integral corresponding to Hq. 
Defining Pf{W, x) as the joint probability that work W is done in the forward process with a 
final value of x for the reaction coordinate, and Pr{W, x) as the joint probability for work W 
and initial value x i^i the reverse process, then one can derive analogously to Eq. flT7|) that 
Pr{—W,x) = e~^^e^'^^P/(iy, x), again using Liouville's theorem. Integrating this identity 
over W gives 

{Six - X[x]))xs = {Six - X[xir)])e~^^^-^^^)x.. (20) 

Since the objective is to get the potential of mean force, the appearance of AF seems to 
be a complication. However, the left hand side of Eq. fl2Ul) does not correspond to the 
potential of mean force in Eq. (TTOl) . since the Hb and Hq ensembles are different. Instead, 



from Eq. (ITSI) . one easily shows that {5{x — x[^]))xb — '§^^ '^'^^^'^^H^ix ~ x[^])) o, and since 
g-/3AF ^ Zb/Za, Eq. ([20]) becomes 

,-/3/(x) = c(<5(x - xk(r)])e-'^[^-*(>^'^-)l>,^ , (21) 



which is the path-integral analog of the result of Paramore et al. 25|. Note that c = -^ 
is a constant, which does not matter for the potential of mean force, so that /(x) can be 
determined from a non-equilibrium process in a similar way as the free energy IS.F . 



IV. REGULARIZATION METHODS 

The invariance of the volume element T>xVp under Hamiltonian dynamics is essential to 
derive the Jarzynski equality in flT6l) . the Crooks fluctuation relation in flT7|) and its exten- 
sion to constrained ensembles. However, strictly speaking the phase space volume element 
is infinite here. The infinite volume of the phase space volume element is reminiscent of ul- 
traviolet divergences in classical field theories that arise from an infinite dimensional phase 
space. For instance, the average kinetic energy of a classical elastic string is equal to the 
number of degrees of freedom times 1/(2/3), but since the number of degrees of freedom of the 
string is infinite, the average kinetic energy diverges. The ultraviolet divergences present 
difficulties in the direct application of the results of Scholl-Paschinger and Dellago[2l[ to 
the elastic string, because in the definition of the work (TT5|) . both H[x{t),p{t), Xb] and 
H[x{0),p{0), Xa] are divergent quantities, so W might not be well defined. Other compli- 
cations would arise for dynamics with phase space contraction, which are not considered 
here. 

To assess whether the fiuctuation relations derived in the previous section are meaningful 
for a system with an infinite number of degrees of freedom, the divergences need to be 
regularized such that a finite number M of degrees of freedom results. There are two general 
approaches to regularizing the path integral: the first is a Fourier regularization, in which a 
wave vector cut-off in Fourier space is introduced, while the other is a bead regularization, 
which consists of discretizing the points of the elastic string by replacing the continuous 
imaginary time s by a finely spaced lattice of imaginary time points. In order to establish 
a proper theory, the regularized quantities such as the free energy and the distribution of 
fictitious work must converge to a finite limit as the M — >■ cxd, corresponding to the imaginary 
time lattice spacing going to zero or the cut-off to infinity, respectively. Such regularization 
is also required to obtain feasible computational methods, and for these the nature of the 
convergence (of the free energy, the work distribution etc.) to a finite result is important 
for the efficiency of the method. 

While taking a finite value for M solves the infinite phase space volume problem for 
the partition function using either of the two regularization procedures, it is a separate 
question whether the distribution of work values is well defined in the limit M -^ oo. The 
work distribution is used in the Jarzynski and Crooks relations, and is central to the non- 
equilibrium methods. For any finite M, this distribution will be well defined and yield 
information on the finite-M free energy difference AFm- While Uthm^oo'^Fm is equal 
to the true quantum free energy difference AF, it must be noted that there is currently no 
general method to show that P{W) is well behaved as M — >■ cxo. In Sec. |Vl it is demonstrated 
numerically that the work distribution converges for an asymmetric double-well potential. In 



the companion paper [18|, the convergence of the work distribution will be shown analytically 



for the specific case of a particle in a harmonic well of changing strength, which can be solved 

exactly. 

A. Fourier regularization 

The central idea of the Fourier space regularization approach is to restrict the number of 
Fourier components x^. and pk of the continuous fields x{u) and p{u) to be finite, where 

Xk = {Tx)k (22a) 

Pk = (^p)fc, (22b) 

and 



Jo 



'0 

Note that because x and p are periodic with period 1, k only takes integer values. Further- 
more, since x and p are real fields, their Fourier modes satisfy x_fc = xl and p^k = Pl- In 
the Fourier representation, the Hamiltonian flT^ becomes 

H{±, p. A) = f; (^ + ^mcollikl') + Vi±, A), (23) 

fc=— oo 

where x and p are the collections of all Xk and pk, respectively, and the dispersion relation 
is 

ujk = 27rkJ- = -— 24 

V m Tip 

and 

V{5^,X)=J^{V{J^-%X))k=o. (25) 

Note that for A; 7^ modes k and —k are degenerate. Given the Taylor series for the 
potential, V{x, X) = ax + h\x + cxx"^ + dxx^ + . . . , from Eq. fl25|) one gets 

V(x, A) = aA + ^A^o + CA ^ \Xk\^ 

k = — OD 

00 00 

+ ^^ 5Z 5Z ^'fci^fca^fci+fca + • • • (26) 

fci= — 00 fc2 = — 00 

Because of the periodic boundary conditions imposed on the fields x{u) and p{u), the 
infinite volume in phase space is now countable. We can thus regularize the theory by 
imposing a cut-off kc on the values of k. With a finite cut-off, the application of Liouville's 
theorem to derive Eqs. (TT6!) and (TT7|) poses no problem since the Hamiltonian flow involving 
M = 1 + 2kc degrees of freedom preserves phase-space volume for any finite M, and hence 
also in the limit M — > 00. Note that the limit M — ;> 1 corresponds to the classical limit, as 
can be seen by taking only the k = term in Eq. fl23j) . 

Using this Fourier space regularization, the free energy converges in the limit M — >■ 
cx) as 0{M~^) when used in the above straightforward form 26l|. Using so-called partial 



averaging techniques [2 7i], this can be turned into a 0{M~'^) convergence [26|]. However, the 
next regularization approach, based on replacing the string by a set of beads connected by 
springs, is a much simpler and more general way to get an 0{M~'^) convergence (or higher). 

8 



B. Bead regularization 

In contrast with the Fourier regularization, in the bead regularization procedure, the 
closed path x{u) is represented by a lattice of M points. The points u = n/M of the lattice 
are called beads, and their positions and momenta are denoted by x = {xi, . . . ,xm} and 
P = {Piy ■ ■ ■ ,Pm}, respectively. Using these variables, the partition function Z{X) in Eq. fITT]) 
can be approximated byfl7| 



Zm{X) 



mM 



M/2 



/"dxdp e-^^^^^'^'P'^), 



where 



M 



H^ 



M 



;x,p>^) = Y1 



n=l 



P^ 1 

2/i' 2 ^ " 



Xn+lf + J^V{Xn,\) 



(27) 



(28) 



and Xm+1 = Xi, i.e., the x„ form a "ring polymer." In Eqs. fl27|) and fl28|) . fi' is the mass 
associated with each bead. To ensure that the ring polymer approaches the elastic string 
limit as M — » oo, the mass of each of the M beads has to be fi' = fJ^/M, but for numerical 
applications at finite M, n' is a free parameter. 

This regularization scheme is the basis of the frequently used path integral molecular 
dynamics (PIMD) method IT], |28|, |29|, |30| for computing canonical equilibrium averages for 
quantum systems. 

To derive the bead regularization of the partition function in Eq. ( 1271) and the fictitious 
Hamiltonian in Eq. fl25]) . one usually starts from the Trotter formula for the Boltzmann 
operator, i.e. 3l|, 1321 



-f3H 



-I3V/M -I3f/M\ 



M 



lim [e '^' '"' e , 



(29) 



However, it is hard to see from Eq. fl2I?l) why the convergence of the free energy for large 
M would behave (9(M~^), as mentioned above. The convergence properties are easier to 
determine starting from the symmetric version 



-I3H/M 



-PV/i2M)-PT/M-l3V/(2R4) 



+ 0{M- 



(30) 



which follows from the Baker-Campbell-Hausdorff formula[33|. Equation fl5D]) allows the 
Boltzmann operator to be expressed as 



-PH 



g-/3y/{2Af)g-/3T/Afg-/3y/(2M) 



M 



+ 0(M" 



(31) 



When taking the trace of the Boltzmann operator to obtain the partition sum in Eq. ([T]), 
the first term on the right-hand side of Eq. (I3T]) can be rewritten using the cyclic properties 
of the trace in the form of the Trotter formula (!29|) . Thus, for the partition sum, the two 
so-called splitting methods (129!) and (13T|) lead to the same result. Since the latter converges 
as 0(M~^), so does the former. Note that if the trace is not taken, such as for the imaginary 
time propagator or for expectation values of operators which do not commute with x, the two 
splitting methods exhibit different convergence behavior 3J] . In particular, the convergence 
behavior of the Trotter form then becomes (9(M~^)|32l|. 



Taking; the trace to get the partition sum, and the standard technique of inserting closure 



relations |15l. Il6l . Il7l |. one finds 

M/2 



Z^iX) = (^) |dxe-^^"-[i^(---^^^^-^^(-)]. (32) 



Finally, multiplying the right-hand side of Eq. (15^ with 

^) /dP»-^^-- (33) 

one obtains Eqs. (l27j) and fl28l) . and one sees that Zm{X) converges to Z{X) as (9(M~^). 

The bead regularization in Eq. fl28|) is not as different from the string's Fourier represen- 
tation in Eq. fl23|) as it may appear at first sight, as becomes clear when Eq. fl28l) is written 
in the Fourier representation as well. Since the beads are discrete, the Fourier transform is 
also discrete and takes the form 

x. = ^5^e2-''^"/^a;. (34a) 

n=l 

M 

Pfc = 5^e2-'^"/^p„, (34b) 

n=l 

where k runs from to M — 1. Because Xk = XM+k, one can also choose the range of k to 
be from — [M/2J to [M/2J , and this is more convenient since a natural cut-off kc = \_M/2\ 
then arises. The asymmetry between the position and the momentum transformation in 
Eqs. (I34aj) and (134bp is necessary to have the Fourier transformation preserve the canonical 
structure, while at the same time letting the definition of Xk in Eqs. f l22al) and fl34ap coincide 
for M ^ oo. Applying this transformation to the fictitious Hamiltonian in Eq. fl28|) gives 

where the only real difference with the fictitious elastic string Hamiltonian (!23ll (with a wave 
vector cut-off kc) lies in the dispersion relation 

2M . tt/c 

instead of Eq. (l24l) . Note that for k ^ 0, k and —k are again degenerate, with the exception 
that for even M, the mode k = —M/2 is identical to the mode k = M/2 and only one of 
these should be included. Naturally, in the limit M — > oo, Eq. (136|) reduces to Eq. (^^ for 
fixed k. 

The potential term V in Eq. (!35l) is given by the same expression fl26l) as for the elastic 
string, with all sums over wave vectors cut-off at kc- Therefore, both for the Fourier and the 
bead regularization, the equations of motion in terms of x„ and p„ take the form 

dxk _pk_ 
dt ~/i" 

-— - = - ulxk - bx5ko - 2cxXk - Sdx ^ XgXk^g - ■■■ (37b) 

q=-kc 

where ii" = n for the Fourier regularization and /x" = M/x' for the bead regularization. 

10 



(37a) 



V. PATH INTEGRAL SIMULATIONS 

As an illustration of calculations of free energy differences in non-trivial quantum systems, 
consider the difference in free energy between a quantum particle confined in a symmetric 
quartic double- well potential Va{x) and a quartic potential with a linear bias Vb{x) (see 
Fig. DP) 

Va{x) = Vo {x^ - x^) (38a) 

Vb{x) = Vo {x* -x'^ + x). (38b) 

The free energy difference for a quantum particle confined by the potentials Va and Vb 
will be computed from the Crooks fluctuation relation using the time-dependent potential 
V{x, X{t)) = Va{x) + \(t)Vox, where we assume A(t) = t/r and r defines the rate at which 
the potential is switched. Note that Xa = A(0) = and A^ = A(r) = 1. The bead 
regularization will be used because of its better convergence properties compared to the 
Fourier regularization. 

For application of the non-equilibrium relations, the procedure consists of drawing A^ ini- 
tial values of the bead positions x = {xi, . . . , xm} and conjugate momenta p = {pi, . . . , pm} 
from the canonical probability density 

p(x, p) = e-^^"("'P'°V^A 

and propagating each phase point (xj, pj) {i = 1 . . . N) forward to time r under the influence 
of the time-dependent potential V^(x, A(t)). 

A. Sampling initial conditions 

The sampling of initial conditions for the dynamics was done using a Monte Carlo pro- 
cedure. For systems that have weak quantum character, i.e., with large k = m/{(3hy, the 
harmonic part 

f//,(x) = -kM ^ {xn - a;„+i)^ 

n=l 

of Hm arising from the kinetic energy operator of the quantum particle, forces the bead 
positions Xn to be near one another. This strong harmonic potential Uh{x.) makes Monte- 
Carlo sampling of the canonical distribution inefficient if the sampling is based only on 
trial moves generated by uniformly chosen random displacements of the bead positions Xn- 
For this reason, it is helpful to use importance sampling based on the probability density 
for a free particle. The procedure consists of dividing the trial configurations generated in 
the Monte-Carlo sampling into two types. The first type consists of randomly displacing 
the centroid or center-of-mass of the ring polymer. Since the potential Uh is constant for 
this type of displacement, only the potential \^(x. A) determines the acceptance probability 
of this type of move. The second way of generating trial configurations consists of drawing 
independent, non-zero Fourier modes {xk} with k = Ito k = [^^y^\ based on the probability 
density 

P(4),4*)) = ^e-'^^-H'^+^n, (39a) 

71 
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FIG. 1: Two model potentials confining a one dimensional quantum particle, for which the free 
energy difference is determined by non equilibrium methods in Sec. |Vl Note that the A values 
given correspond to V{x, X = 0) = Va and V{x, A = 1) = Vg. 



where a^ = Prnul and x^ 



that if M is even, then x^L^ 



~(r) 

and X 



and x)^' are the real and imaginary parts of x^, respectively. Note 



M/2 



is drawn from 



^y^M/2> 



)l 1. -(r)2 



(39b) 



with bk = ^PmLul. The Fourier modes {xk} are then inverted to form the positions of the 
beads x for the trial configuration, which is then accepted with probability min (l, e~^^^), 
where AV is the difference in the potential between the trial and original positions (i.e. 
with no harmonic contribution Uh)- Since this procedure generates trial configurations 
with the same centroid as the original configuration, A^ is typically small and most trial 
configurations are accepted if quantum effects are not too large. Once a statistically inde- 
pendent configuration of the ring polymer has been obtained, conjugate momenta p may be 
drawn from a Maxwell-Boltzmann density to obtain an initial phase point for the system 
[cf. Eq. m]- 



B. Dynamics 

Although the numerical time propagation of the system can be carried out in a multitude 
of ways, symplectic integration methods usually offer superior stability and accuracy. The 
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general Hamiltonian (l28l) is time-dependent and the dynamics is more complicated than that 
of autonomous systems. Nonetheless, it can be shown that phase space volume still satisfies 
Liouville's theorem under the dynamical fiow and that integration schemes of second and 
higher order can be derived for the time-dependent potentials in the same manner as for 
autonomous systems, as long as the time is updated only after the momentum-propagation 
step(35|]. For our test system, we used a second-order Verlet-like integration scheme that 
satisfies this requirement, where first the momenta are propagated for a half-step 5t/2 using 
initial forces, then an update of the positions and the system time is done using the current 
momenta by a full step 5t. The time-dependent force is then computed at this system 
time using the updated positions, and finally a final momentum update of a half-step 5t/2 
is carried out. It is readily established that a trajectory using this scheme is exact for a 
Hamiltonian that differs from Hm by terms of order 5t^ at all times. 



C. Estimating the free energy difference 

After each phase point has been propagated to time r, the non-equilibrium work Wi = 
HMi^i{T), Pi{T), 1) — if (xj(0), Pj(0), 0) is computed. From a set of n initial conditions and 
work values, the free energy difference for the quantum system can then be computed using 
the Jarzynski estimator 



-/5AF = logfif]e-^->). (40) 



=1 



Statistical uncertainties for this estimator may computed using jackknife [36!] or bootstrap[3 



38j methods on the sample. 

As is clear from Eq. fHUl) , the free energy difference computed from the exponential average 
of the work is sensitive to large fluctuations in the value of the work Wi and may converge 
slowly if there are large tails in the work distribution|39i]. The Crooks fluctuation relation, on 
the other hand, is based on flnding the value of the work Wc at which the probability density 
of the work Pf{Wc) is equal to the probability density Pr{—Wc) of the negative of the work 
in the reverse process in which initial conditions are drawn from a canonical distribution 
based on the potential Vb(x). According to the Crooks relation in Eq. fITTj) . the free energy 
difference is then given by AF = Wc- Calculations of the free energy difference based on 
the Crooks fluctuation relation therefore require constructing the probability densities of 
the work in the forward and reverse directions. The traditional approach of approximating 
such densities is to use histograms. While the use of histograms is parameter free for 
systems in which the variable is conflned to discrete values, the more typical situation 
concerns probability densities of continuous variables, which requires the speciflcation of 
a bin size for representing the density. Thus, in contrast to the calculation of the free 
energy difference using the Jarzynski fluctuation relation, the use of a histogram approach 
to approximating probability densities leads to an undesirable parameter dependence to the 
free energy differences computed in the Crooks fluctuation approach. 

In fact it is possible to reconstruct probability densities of continuous variables without 



resorting to parameter-dependent histogram methods |40[|. The idea is to expand the em- 
pirical cumulative distribution function (ECDF) obtained after sorting the data in a series 
of complete orthogonal polynomials. From the mathematical properties of the cumulative 
distribution function, the number of terms required in the expansion of the ECDF can be 
determined without user intervention by application of the Kolmogorov or Kuiper's test 38| . 
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From this expansion, an analytical form for the probability density can be obtained by dif- 
ferentiation. This approach can be applied to the probability densities of the work for both 
the forward and reverse processes. The free energy is then obtained numerically by finding 
the value of the work where the analytical probability densities are equal. 

In principle, any orthogonal set of polynomials can be utilized for the expansion of the 



ECDF. In Ref. (l40i). a Fourier series for the difference between the ECDF and a linear 

function was used to construct an analytical approximation to the ECDF. Although such 

an approach successfully produces a fairly smooth approximation of the ECDF, the use 

of oscillatory transcendental functions introduces high frequency oscillations in the smooth 

approximation that may lead to systematic errors when finding points of overlap of densities, 

particularly when the densities intersect in regions where there are prolonged tails. As 

an alternative, we consider using a high-order polynomial expansion via the Chebyshev 

polynomials. 

To illustrate the method, consider a series of n sorted work values {wi} where Wi < ti^j+i. 

An unbiased estimator of the distribution function of the work P{w) is the ECDF defined 

in the range [wi , ti;„] as 

_ ^ 

P{w) = — ioT Wi < w < Wj+i. (41) 

n 

To allow an expansion in Chebyshev polynomials which are defined on the interval [—1,1], 

the ECDF is mapped to that domain using 

_ 2W-Wi- Wn ,.^. 

w = (42) 

Wn -Wi 

F{w) = P{w). (43) 

The new ECDF F can then be expanded in terms of Chebyshev polynomials of the first 
kind T„ to yield the approximation 

Piw) ^ PUw) = FUw) = - + - E djT.iw), 

where the coefficients dj are given by 



-1 V 1 — w"^ 



"3 



The analytical approximation to the probability density Pm{w) in terms of m Chebyshev 
polynomials is then given by 



4 



Pm{w) = — Vjrfjf/,_i(w), 

7T{Wn-Wi) ^ 



where f/„ are Chebyshev polynomials of the second kind|4ll]- The expansion coefficients d 



^3 



can be evaluated analytically using the form of the ECDF in fl4T|) . and are given by 



1 " 
do = —y axcco^iwi 
n ^-^ 



1 " 1 I 



n -^ — ' ? 
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In practice, one hopes that the number of polynomials m required in the expansion of the 
ECDF is modest so that a smooth approximation is obtained. What number of polynomials 
is appropriate can be estimated using either the Kolmogorov|38l. |40| or Kuiper's test(38|, 



which determine how likely it is that the difference between the ECDF P{w) and its an- 
alytical approximation Pm{w) is due to random variations. The tests take the maximum 
variation between P and P^ over the sampled points and return a probability Qm that 
the difference between the two cumulative distribution functions is due to chance. A small 
value of Qm indicates that the difference between the cumulative distribution functions is 
significant, so that the quality of the expansion Pm is insufficient to represent the data. 
One therefore carries out a process of progressively increasing the number of polynomials m 
and evaluating Pm as well as Qm until the value of Qm is larger than some threshold, say 
Qc = 0.5. 

Once analytical approximations to the probability densities of the work in the forward 
and reverse directions have been obtained, the intersection point Wc is readily evaluated by 
numerically searching for a solution of PfiWc) = Pr{—Wc) using the Brent method[38|. The 
Brent method requires that the solution be bracketed on an interval {wmin,Wmax)- Since 
the free energy difference between ensembles are single-valued, there is only one point of 
intersection of the probability densities in the forward and reverse directions. The interval 
end-points Wmin and Wmax can therefore be taken to be near the extremal values of the Wi 
values found in the simulations. 

Finally, statistical uncertainties can also be computed for the free energy difference by 
repeating the calculation of the intersection point for a series of jackknifed samples of the 
data and using the variance of the free energy differences over the jackknife samples. 

D. Simulation results 

To illustrate this approach, numerical tests were carried out to calculate the quantum 
free energy difference between systems with A = (symmetric quartic double-well potential) 
and A = 1 (biased quartic potential) with parameters (3 = l,h = l,m = l and Vq = 5, while 
the mass /i' per bead was also set to 1. The time step for the propagation was 6t = 0.001 
such that the fluctuations in the total energy at constant A = relative to the fluctuations 
in potential energy were less than 1%. The calculations were carried out at two different 
switching rates, r = 0.5 (fast switching) and r = 100 (slow switching). 

In Fig. [21 the free energy difference AFm calculated using the Crooks fluctuation relation 
for fast- switching process is shown as a function of the number of beads M in the ring poly- 
mer regularization of the path integral. The work value Wc at which Pf{Wc) = Pr{—Wc) 
was based on the analytical approximation of the empirical cumulative distribution func- 
tions formed out of 1 x 10^ independent realizations of the non-equilibrium process in the 
forward and reverse directions. The inset graph in this figure clearly suggests that AFm 
converges as 1/M^, to a final value of AF = —2.35. This value agrees with the quantum free 
energy difference found by evaluating the partition sum using the numerically determined 
eigenvalues of the Hamiltonians H{0) and H{1)- The free energy difference is not just due 
to the difference in zero-point energy, which is —2.53, but to the higher energy levels as well. 
Note also that the quantum contributions to the free energy difference lead to a quantum 
free energy difference that is roughly 25% higher than the classical value of —2.95. 

In Fig. [31 the first and second cumulants of the work done (i.e., the average and the 
variance of the work) in the fast reverse process are plotted against M. As the fits in the 
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FIG. 2: The difference in free energy for the quartic system as a function of the level of discretization 
of the path integral for the fast-switching non-equilibrium process (r = 0.5). The solid line is drawn 
as a guide to the eyes, while the dotted line in the inset is a fit to a + 6M~^, as predicted by the 
theory, which agrees within errorbars. 

figures show, the asymptotic convergence of both cumulants is consistent with a 0{M^'^) 
behavior. The results presented in this figure suggest, in general, that the work distributions 
converge in the infinite M limit as (9(M~^), which can be rationalized from the consideration 
of a harmonic oscillator system (see the companion paper [18|). 

The probability densities of the work and negative work in the forward and reverse di- 
rections are shown in Fig. [Has histograms and expanded analytical forms. The convergence 
properties of the analytical expansions of the empirical cumulative distribution function were 
determined by the asymptotic Kuiper's test with a threshold value Qc of 0.5. Typically, 
m = 13 terms were required to reach convergence with this choice of Qc- The probability 
densities shown in Fig. H] were estimated from 1 x 10^ values of the work and negative work 
for the non-equilibrium switching process. The data shown are for a M = 9 bead discretiza- 
tion of the path integral. Note that even though the shape of the probability densities is 
sensitive to the switching rate, the intersection point of the forward and reverse densities is 
independent of the rate and equal to the free energy difference. It is apparent from the detail 
of the work densities near their point of intersection shown in the inset of the panels that it 
is difficult to arrive at an estimate of the free energy difference based on noisy histograms 
of the work. In contrast, the analytical forms of the probability densities lead to smooth 
curves and unambiguous points of intersection. 
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FIG. 3: The first and second cuniulants of the work distribution for the reverse process for the 
quartic system as a function of the level of discretization of the path integral for the fast-switching 
non-equilibrium process (r = 0.5). The solid lines are a fit to a + bM^^. 

VI. CONCLUSIONS 

Non-equilibrium methods for the calculation of free energy differences in quantum systems 
in the context of the path integral representation of the canonical partition function have 
been presented. Instead of using the real quantum dynamics of the system, the path integral 
representation allows a fictitious path to be defined for which the Jarzynski and Crooks 
relations are valid. By evolving the ring polymer in the path integral representation under 
fictitious dynamics, the difficulties associated with the complexity of the full evolution of 
a quantum system are avoided. From a computational perspective, avoiding true quantum 
dynamics is a great advantage, although the actual efficiency as compared to other methods 
to compute free energy differences will depend greatly on implementation details and on 
using additional techniques such as importance sampling. 

From a formal point of view, the path integral approach exploits the well-known isomor- 
phism between a quantum system and a classical system of higher dimensionality described 
by a field theory. While the dynamical evolution of the classical field is well-defined, the 
isomorphic classical system exhibits the typical divergent behavior of classical field theories 
on short length scales, due to the infinite number of degrees of freedom. It was demonstrated 
that the equations of motion and field variables can be regularized using either a wave- vector 
cut-off of the Fourier modes of the fields x{u) and p{u), or a real space discretization of the 
ring polymer representing the quantum particle, yielding a finite number M of degrees of 
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FIG. 4: (Color online) The forward and backward work distributions in histogram and analytical 
expanded form based on 1 x 10^ realizations of the non-equilibrium process for the quartic system, 
using M = 19. The top panel corresponds to results for the fast-switching process with r = 0.5 
while the lower panel contains the results for a slow-switching process with r = 100. Noisy curves 
correspond to the histogram form and the smooth curves to the analytical expanded form of the 
work distributions. The insets show the detail in the vicinity of the crossing point. 

freedom. 

A general numerical procedure for calculating the free energy difference between non- 
trivial quantum systems was elaborated using a particle confined in a quartic potential as a 
test model. General issues pertaining to sampling the initial phase points, performing the 
dynamics of non-autonomous systems, and estimating free energy differences using the expo- 
nential average of the work ( Jarzynski method) and the crossing method (Crooks approach) 
were discussed. A parameter-free method for calculating the free energy difference using the 
crossing of the forward and reverse work distributions was introduced. The parameter-free 
method is based on expanding the empirical cumulative distribution function in orthogonal 
Chebyshev polynomials and is controlled by a rigorous statistical convergence test. From 
these expansions, analytical functions of the estimated forward and reverse work distribu- 
tions are then obtained from which precise values of the crossing point and hence the free 
energy can be extracted. 

The expansion procedure utilized here is quite general, and can be used to construct 
analytical estimates of any probability density of a single continuous variable. The approach 
should be particularly useful whenever the specific value of a probability density if desired, 
such as in the calculation of the value of the radial distribution function at contact in a 
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system of hard spheres. The method could also be of use in construction of the free energy 
as a function of a reaction coordinate or other parameter, where the dispensing of parameter- 
dependent histograms is desired. 

The numerical results indicated that for the quartic potential the work distribution con- 
verges as the regularization parameter M gets large, and hence so does the free energy 
estimate obtained from the crossing method. Although the general convergence of the free 
energy and work distribution utilized in the crossing method must be demonstrated on 
a case-by-case basis, the following paper shows that for the harmonic oscillator the work 
distribution and free energy converges rigorously but that the nature of the convergence de- 



pends on the regularization used[18|]. In particular, it is shown that the bead regularization 
converges faster than the Fourier regularization. 

While the focus here has been on deterministic methods with Hamiltonian dynamics in 
non-equilibrium statistical mechanics, it is worth mentioning that an alternative to the use 
of Hamiltonian dynamics is to evaluate the path integral via the path integral Monte Carlo 
(PIMC) method. Both the Jarzynski and Crooks relations hold in this context provided that 
the process is Markovian and time reversiblejg, |7|, |8|. Furthermore, for a given regularization 
at finite M , one can in principle also use thermostatted deterministic dynamics, for which the 
Jarzynski and Crooks relations also hold 2l| . The convergence properties of these alternative 



non-equilibrium processes as M ^ oo will be assessed in future work. 

Finally, it should also be noted that the dynamics generated both through deterministic 
and stochastic evolution is completely artificial and generally has little to do with the real 
time quantum dynamics of the system (except when hj3 -^ 0). Still, the computation of 
the free energy difference through the dynamical non-equilibrium procedure described above 
yields exact quantum results in the limit M — > oo. Moreover, the use of either deterministic 
or stochastic evolution bypasses the problem of computing the real time quantum dynamics 
and illustrates the power of path integral methods in practical applications. 
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